. ii , o o h '9 \* .1 1*0^0 if 92 9 1 

JG13B@c'dPCT/PT0 1 1 FEB 2002 



ELECTROMAGNETIC SIMULATION ALGORITHM, IN PARTICULAR FOR 
THE PERFORMANCE OF AN ANTENNA 

The present invention pertains to an electromagnetic 
5 simulation algorithm, in particular for the performance 
of an antenna, which makes it possible to compute the 
electromagnetic wave scattered by a conductor in a 
monof requency situation. It applies in particular to 
the simulation tools used during the design of 

10 reception or transmission antennas such as cellphone 
antennas, anti-collision radar antennas, those of 
electronic counter measures (ECM) systems, of 
monitoring or tracking radars, or satellite antennas. 
The invention also applies to the computation of radar 

15 cross sections (RCS) of objects whose geometrical 
properties are known. 

Antenna simulations are used to limit the number of 
mock-ups and prototypes during the design of said 
20 antennas. These simulations make it possible in 
particular to compute the far-field radiation pattern 
of the antennas and to adapt the antennas in 
transmission or in reception, in the present or 
otherwise of a surrounding structure. As input data 

2 5 they use a mesh of the antenna whose performance one 

wishes to evaluate, as well as the characteristics of 
the electromagnetic excitation to which it is subject. 
The invention is not limited to the simulations of 
antennas. It applies also for example to the 
30 computations of RCSs of targets. The application of the 
invention within the antenna simulations, during 
reception, will now be described by way of 
illustration . 

3 5 Two main methods can be distinguished within the 

simulations commonly employed. A first method is based 
on computation by finite differences, also known as the 
volume finite element method. According to this method, 
a mesh of a volume surrounding the antenna is used. A 



drawback of this method is that the mesh is necessarily 
bounded, whereas one is interested in the radiation 
pattern at infinity. A compromise must then be made 
between the dimension of the meshed volume, that is to 
5 say the accuracy of computation, and the computation 
time. To alleviate this drawback, a second method is 
used, based on integral equations within the frequency 
domain. According to this method, a surface mesh of the 
antenna only is used. The radiation pattern at infinity 
10 is computed directly from electric and magnetic 
currents on the surface of the antenna. 

Certain known techniques using integral equations 
computed the electric and magnetic currents (from which 
15 the field pattern radiated at infinity is deduced) by 
virtue of a factorization of a matrix. This matrix is 
known as the interaction matrix or else the impedance 
matrix. This factorization allows direct computation, 
that is to say noniterative computation, of the surface 
20 currents. A drawback of these techniques is that the 
computation time is long. If N denotes the number of 
points involved in meshing the antenna (also termed 
surface triangulation) , the computation time according 
to these techniques varies as N 3 . Now, the number of 
25 points N is itself related to the wavelength (and 
consequently the frequency) of the wave radiated by the 
antenna. Let us assume that a simulation is carried out 
at 10 GHz, by using N points in the mesh of the 
antenna, and that the computation time is x. To 
30 transpose this simulation to 20 GHz, it will be 
necessary to use a mesh comprising 4xN points, and this 
will represent a computation time of the order of 4 3 xx. 
A computation time problem arises also when one seeks 
to simulate complex antenna geometries, such as small- 
35 arrays. This makes these techniques unusable in 
particular in design tools which require a restricted 
computation time so as to allow the designers to carry 
out several tests . 
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Other known techniques using integral equations make it 
possible to reduce the computation time by virtue of an 
iterative solution method. If IT denotes the number of 
iterations, the computation time according to these 
techniques varies as ItxN 2 . One problem with these 
techniques is that nothing guarantees the convergence 
of the computations. Stated otherwise, antenna shapes 
exist for which the radiated field pattern cannot be 
computed with these techniques. 

An aim of the invention is to alleviate the aforesaid 
drawbacks, and in particular to restrict the 
computation times. 

To this end, the invention relates to an algorithm for 
simulating the performance of an antenna, based on 
iteratively solving a system of integral equations 
comprising a precondi tioner . This precondi tioner arises 
in particular from adapting Calderon's formulae to the 
boundary integral equations of electromagnetism. In 
particular, in the case of a completely metal antenna, 
a preconditioner is proposed for the equation known as 
the "Electric Field Integral Equation" (EFIE) . Use is 
also made of an original representation of the residual 
of the computations during each iteration. This 
representation, as well as a projection and a 
composition, are involved in the expression of said 
preconditioner . 

The invention has the following main advantages: 

• it converges rapidly; 

• it makes it possible to simulate arbitrary geometries 
and any excitations; 

• its conditioning is independent of the fineness of 
the mesh; 

• it accommodates algorithms based on the computation 
of an impedance matrix, by reusing said impedance 
matrix; 
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• it makes it possible to deal with antennas 
comprising, in addition to metal, dielectric 
materials . 

5 Other characteristics and advantages of the invention 
will become more clearly apparent in the description 
which follows and in the appended figures which 
represent : 

- figure 1, a mesh of an antenna; 

10 - figure 2, a sectional view of the mesh of 

figure 1 ; 

- figure 3, a detail of the mesh of figure 1, 
in which a vector field is represented; 

- figure 4, a functional diagram of an 
15 iterative algorithm; 

- figure 5, a trianglewise constant basis 
function over a mesh; 

- figure 6, a trianglewise affine and 
continuous basis function over a mesh; 

20 - figures 7 and 8, two illustrations of the 

performance of the algorithm according to the invention 
as compared with known techniques. 

Reference will now be made to figures 1 and 2 which 

2 5 represent an exemplary antenna shape for which one 

seeks to determine the scattered field when the antenna 
is illuminated by an incident wave. Stated otherwise, 
one seeks to simulate an antenna during reception. The 
antenna taken in this example is a spherical cavity 
30 with aperture half -angle 71/ 4. The inside radius is 7/8, 
the outside radius is 9/8 (arbitrary unit of length) . 
The surface of the antenna is denoted ^ . This surface 
r is meshed with triangles. The surface T is assumed 
to be that of a perfect conductor (the antenna) £2- 

3 5 immersed in a vacuum £2+ . 

In this example, the incident electromagnetic wave 
which illuminates the antenna is a monof requency plane 
wave. This incident electromagnetic wave, of known 
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wavenumber k inc , is represented by two vector fields 
denoted E ,nc and H mc corresponding respectively to the 
electric field and to the electric field and to the 
magnetic field. Of course, the invention does not apply 
5 only to plane waves. It is possible to substitute this 
incident wave with the field emitted by a radiating 
dipole for example (transmitter mode of the antenna). 

One seeks to determine the electromagnetic wave 
10 scattered by the antenna at infinity, that is to say 
the far-field radiation pattern. This scattered 
electromagnetic wave is represented by two vector 
fields denoted E mf and H d,f corresponding respectively to 
the electric field and to the magnetic field. 

15 

The electromagnetic field radiated at any point in 
space can be computed from the field of current flowing 
around the surface of said antenna, also referred to as 
electric and magnetic surface currents. The asymptotic 

20 expression at infinity for the electromagnetic field 
radiated is the radiation pattern which one seeks to 
determine. This computation, well known in 
electromagnet ism, is recalled in the document "Integral 
Equation Methods in Scattering Theory" by 

25 D. Colton and R. Kress - John Wiley & Sons, New York, 
1983 . 

The electromagnetic field satisfies Maxwell's equations 
within the vacuum £2 + which may be written: 
cui-i(e)= icofaH (1) 

c"u7i(h)= -icoeE (2) 

The terms used in these equations represent: 

• curl the curl operator 

• E = E"' c + E"' 1 the total electric field in complex 
notation : 

3 5 • H = H ,nc + H''' / the total magnetic field in complex 
notation : 

• 1=7^1; 
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• 0) the angular frequency of the electromagnetic wave; 

• 11 the magnetic permeability ; 

• e the electric permitivity. 

It is recalled that the wavenumber is related to the 
angular frequency simply by the following relation: 

k = Giy/ile (3) 

The electric field of the scattered electromagnetic 
wave E d,J is expressed at any point B of £2+ on the basis 
of the surface currents denoted u by the following 
relations : 



The terms of the relations (4) and (5) represent: 

• div A ( u ) the divergence taken at a point A of the 
vector field u ; 

• grad B the gradient taken at a point B; 

• dS A a differential surface element ; 

• G k the standard Green's function; 

• AB the distance between the points A and B; 

• k the norm of the wavenumber defined by relation (3) . 

A person skilled in the art will be able to investigate 
other technical elements relating to this computation 
in the abovementioned document insofar as the latter 
forms an integral part of the description. 

The computation of the surface currents, that is to say 
of the vector field u , is determined from the following 
variational equation: 
Vv m(u,v)=l(v) (6) 

in which 



E** (B) = kjG k (A,B) u(A) dS A + ^grad, 




(4) 



(5) 
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(7) 



(8) 

This variational equation (6) is known as the boundary 
integral equation of electromagnetism, or else the 
"Electric Field Integral Equation" (EFIE) . 

In order to solve this variational equation (6) in a 
numerical antenna simulation, we must approximate the 
solution u in a space of finite dimension, the so- 
called discretization space. This space contains vector 
fields which represent surface currents . The dimension 
this space is the number of components required to 
fully describe said vector field, at every point of the 
surface T . This surface being meshed, the number of 
components serving to describe u will depend in 
particular on the number of points of the mesh N, as 
well as the nature of the basis functions serving to 
described the vector field (for example, linear 
functions or functions of degree 2). In the subsequent 
description, we shall by way of illustration take the 
Raviart-Thomas space of lowest degree over the mesh of 
T. This discretization space, which represents surface 
currents, is denoted Dh . The symbol h represents the 
characteristic length of the mesh, or else the accuracy 
of the meshing. Specifically, the dimension number of 
this space Dh depends on the number of points N of the 
mesh, which itself depends on the accuracy h of the 
meshing . 

Reference is now made to figure 3 to describe a basis 
of the space Dh . In this basis, the surface current 
field u is represented by the coefficients of a vector 
denoted U. 
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m(Q, v) = k JjG k (A.B) u(A) ■ v(B) dS A dS B 
- £tt r G * ( A - B > div A<u) div B (v) dS A dS B 

and 

l(v) = - jE ,nc (A) - v(A) dS A 
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The basis of the space of surface currents contains 
elements denoted by (p ± where i is an integer index 
associated with an edge of the mesh of the surface T . 
5 These elements are vector fields defined over the mesh 
of the surface T . The vector field (p i , represented by- 
arrows in figure 3 , has a support bounded to two 
triangles Ti and T 2 of the mesh. These triangles Ti and 
T 2 share the edge of index i and of length 1± . This edge 

10 is oriented by a vector of unit norm denoted by y, . Pi 
denotes the vertex of the triangle Ti not contained on 
the edge i ; P 2 denotes the vertex of the triangle Ti not 
contained on the edge i ; P 2 denotes the vertex of the 
triangle T 2 not contained on the edge i . Si and S 2 

15 denote the surface areas of the triangles Ti and T 2 . Let 
z 2 and z 2 be the vectors with unit norm, having a 
direction normal to the surface of the triangles Ti and 
T 2 , and oriented from the interior jQ- to the exterior 
Q.+ . We define the vectors x^ and 3ci, to be the vectors 

20 with unit norms such that the triple (x^,y.,z ) and the 

triple ( x' y.,z ) are right-handed trihedral. We define 
2 1 2 

the vector field (p, for any point A belonging to the 
surface T with the following relations: 



- " AeT 1f cp i (A)=±^-P 1 A with <p, (A) - x} > O ; (9) 

- if A e T 2 . $,(A) = ±^-P 2 A with ^,(A)-x' 2 >0; (10) 

- otherwise, ^,(A)= 6. (11) 



Such a basis of the space of currents is known as the 
usual basis of the Raviart-Thomas space - or else Rao- 
Wilton-Glisson elementary currents. A person skilled in 
30 the art will be able to investigate other technical 
elements in "Electromagnetic scattering by surfaces of 
arbitrary shape" by S.S.M. Rao, D.R. Wilton and A.W. 
Glisson - IEEE Trans. Ant. Prop. AP-30, pp. 409-418, 
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1982 - insofar as this documents forms an integral part 
of the description. 

U± denotes the coefficients of the vector U. The vector 
5 field u(A) used in relation (4) decomposes over the 
usual basis of the Ravi art -Thomas space in the 
following manner: 
G(A)=XU,<p,(A) (12) 

We now transpose relations (6), (7) and (8) into matrix 
10 notation by using the basis ( q>. ) described above. The 

variational equation (6) can be written in the form of 
the following system of linear equations: 

M U = L (13) 

15 

The terms of relation (13) are as follows: 

• M is a known matrix, the so-called interaction matrix 
2 0 or impedance matrix ; 

•Lisa known vector, whose coefficients represent the 
incident wave, that is to say the electromagnetic 
excitation; 

• U is the vector which we seek to determine, whose 
25 coefficients represent the surface currents. 

We define the coefficients M i;) of the interaction matrix 
M and the coefficients Li of the vector L representing 
the incident wave by the following relations: 

30 

M « = k JjG k (A.B)(v,(A)-$ J (B))dS A dS B 

Ber 

-l A JjG k (A.B)div A (ep j )div B (^ J )dS A dS B (14) 



L,= j(E*~(A).cp i (A))dS A 



(15) 



A person skilled in the art will be able to investigate 
other technical elements relating to the computation of 
these coefficients in "Approximation par elements finis 
de surface de problemes de divergence des ondes 
5 electromagnetiques [Finite surface element 

approximation of electromagnetic wave diffraction 
problems]" by A. Bendali - Thesis of the University of 
Paris VI, 1984 - insofar as this document forms an 
integral part of the description. 

10 

Reference is now made to figure 4 wherein is 
illustrated an iterative algorithm for solving the 
system of linear equations (13) based on the conjugate 
gradient technique. It should be noted that the 

15 preconditioning technique which we illustrate in the 
case of the conjugate gradient algorithm applies 
equally well to other iterative algorithms. Mention may 
be made in particular of the Generalized Minimum 
Residual (GMRES) algorithm, Bi Conjugate gradient 

20 (BiCG) algorithm, Quasi-Minimal Residual (QMR) 
algorithm and the BiConjugate Gradient Stabilized (Bi- 
CGSTAB) algorithm. A person skilled in the art will be 
able to investigate other technical elements regarding 
iterative methods in "Templates for the Solution of 

25 Linear Systems : Building Blocks for Iterative Methods" 
by R. Barrett, M. Berry, T.F. Chan, J. Demmel , J. 
Donate, J. Dongarra , V. Elijkhout, R. Pozi, C. Romine 
and H. Van der Vorst - SIAM (1994), Philadelphia, PA - 
insofar as this document forms an integral part of the 

3 0 description. 

Figure 4 illustrates an algorithm 40 which takes as 
input the matrix M and vector L of equation (13) and 
gives as output the vector U. A person skilled in the 
35 art will be able to investigate other technical 
elements relating to the solving of systems of linear 
equations in "Iterative Methods for Linear and 
Nonlinear Equations" by C.T. Kelley - SIAM Frontiers in 
Applied Mathematics, Philadelphia, 1995 - insofar as 



- 11 - 

this document forms an integral part of the 
description . 

A first initialization step 41 makes it possible to 
initialize four series of vectors denoted U[n], R[n], 
S [n] and P[n] where N is an integer. The first series 
U[n] is an approximate solution which converges to the 
sought-after solution U. The second series R[n], known 
as the residual, converges to the zero vector. The 
third series S[n], dubbed the preconditioned residual, 
also converges to the zero vector. The last series P[n] 
is known as the search direction. The first values of 
these series are defined by the following relations: 

U[0]=0 (16) 

R[0] = L - M U[0] (17) 

S[0]=ZR[0] (18) 

P[0]=S[O] (19) 



Relation (16) is used by default when no approximate 
solution is known. A variant of this relation consists 
in taking U[0] to be the result of a surface current 
2 0 computation carried out for one and the same antenna 
but at another frequency. 



The term Z in relation (18) is a matrix. This matrix is 
a preconditioner for the matrix M according to the 
2 5 invention. It is recalled that a preconditioner of M is 
a matrix approximating the inverse of M. According to 
one variant, the preconditioner Z can be replaced by 
the identity. Stated otherwise, S[0] can be initialized 
with R[0] . 

30 

A second iteration step 42 makes it possible to compute 
the values of the aforesaid series at a rank n+1 from 
the terms of rank n. This iteration step uses the 
following relations: 
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U[n+l]=U[n] + aP[n] (20) 
R[n + 1] = R[n] - a M Pjn] (2i } 

S[n + l]=ZR[n + l] (22) 

with 

(R[n],S[nj) 

(MP[n],P[n| (24) 

where (,) represents the complex scalar product in 
5 matrix notation. 

A last step 43 carries out a convergence test. This 
test can be expressed for example by the following 
inequality : 

10 

I R f n J II 

where r\ is a predetermined threshold. Stated otherwise, 
the normalized preconditioned residual is compared 
against the a predetermined threshold r\ . When the 
15 inequality (25) holds, the computation is halted and we 
take U=U[n] as solution. In the converse case, the 
value of n is incremented and we return to step 42. 

Of course, it is possible to use another convergence 
20 test to stop the iterations. It is for example possible 
to replace the inequality (25) by the following 
inequality : 

II sfnJ II 

25 

The conjugate gradient algorithm in the contemporary 
techniques is used without a precondi tioner , that is to 
say with the matrix Z equal to the identity. 
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The invention consists in using a preconditioner based 
on a generalization of a formula of Calderon. A person 
skilled in the art will be able to investigate 
technical elements regarding this formula of Calderon 
5 in "Mathematical Methods in Elec tromagnetism, Linear 
Theory and Applications" by M. Cessenat - World 
Scientific Publishing Co., page 89, 1996 - insofar as 
this document forms an integral part of the 
description. This formula may be written: 

10 

+ J <RJJ 0^= % I (27) 



In this relation (27), 94,J and are three operators 
over the fields of tangent vectors, that is to say 
15 mappings which associate with a field of tangent 
vectors another field of tangent vectors, and I is the 
identity mapping. J is the operator of scalar product 
with the normal to the surface of the antenna, and % 

is the operator associated with the interaction matrix 
20 for which we seek a preconditioner. The operators % j 
and are defined formally by the following relations: 

(5Mu)(B) = [k jG k (A,B)u(A)dS A 

L Aer 

+ ^grad B QjG k (A,B) div A (u) dS A JJ (28) 

(;Q)(B) = Q(B)az(B) (29) 



(<£G)(B) = Jgrad B (G k (A,B))^u(A)dS A (30) 

Aer 

25 In relations (28), (29) and (30), J u , <^ u and u 

are tangent vector fields. The index t in relation (28) 
represents the tangential component of the vector 
between square brackets. B is a point of the surface T. 
z(B) is a unit vector, normal to the surface T at B and 

3 0 oriented outward. 
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The Applicant has found that the operator J ^ J ^ 
being compact, that is to say negligible, the operator 
4 3 94, J is approximately an inverse of the operator 
% . Next, the precondi tioner being defined to within a 

5 multiplicative constant, it is possible to eliminate 
constant 4. Knowing that j* = - j t the preconditioner 
according to the invention is defined from J* J4 J. It is 
recalled that J* is the operator adjoint to J, that is 
to say the operator satisfying the following relation: 

10 

jj*u v = (ujv (31) 

The preconditioner according to the invention is a 
matrix formulation of the operator y* M J. It is in this 

matrix formulation that the advantage of using the 
operator j* % J rather than the operator J <M J 
appears. Specifically, the matrix formulation Z of the 
operator J* 94. J is a symmetric matrix, this being 
essential for iterative algorithms, such as the 
conjugate gradient. 

An exemplary preconditioner according to the invention 
which makes it possible to speed up the convergence of 
the algorithm and also to render this algorithm more 
2 5 stable (we always converge to the solution regardless 
of the initial conditions) will now be described. This 
preconditioner is adapted to electromagnetic problems 
and exploits the structure of the problem to be solved. 

30 We firstly define spaces and their assocated bases 
which will serve subsequently in the description. We 
have already defined the basis ( (p. ) . The space spanned 

by this basis, Dh, is a space of tangent vector fields. 
Stated otherwise, Dh is a set of functions which with 
35 any point of T associate a vector tangent to the meshed 
surface T . We also define two function spaces Ch and 
Sh. The functions of these spaces associate a scalar 
value with each point of T . 



15 



20 
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Ch is the set of trianglewise constant functions, that 
is to say those whose value is constant for any point 
belonging to a given triangle. The basis of Ch which we 
use subsequently is the set of functions denoted Wl, 
which equal 1 on the triangle of index i and 0 outside. 
Reference is made to figure 5 in which a basis function 
W ± is represented. A plain, regular mesh 50 has been 
represented for the sake of clarity. The values of the 
functi on are represented on an axis 51, 

perpendicular to the mesh 50. The function T n 
represented in this figure equals 1 on the triangle i 
and 0 outside. 

Sh is the set of trianglewise affine continuous 
functions. These functions have a constant gradient 
over any given triangle. The basis of Sh which we use 
subsequently is the set of functions denoted 0, , which 
equal 1 at the node of index i and 0 on the other 
nodes. Reference is made to figure 6 in which a basis 
function 0, is represented. A plain, regular mesh 60 
has been represented for the sake of clarity. The 
values of the function 0, are represented on an axis 
61, perpendicular to the mesh 60. The function 0 
represented in this figure trianglewise is affine and 
equals 1 at node i. This function has a support (non- 
zero values) bounded at the triangles 62, 63, 64, 65, 
66, 67 which have a vertex coinciding with node i. 

A person skilled in the art will be able to investigate 
further technical elements in "Handbook of Numerical 
Analysis Vol. II, Finite Elements Methods (Part 1)" by 
P. G. Ciarlet - Ed. J.L. Lions, North-Holland, 1991 - 
insofar as this document forms an integral part of the 
description . 

We shall now define matrices which will serve in the 
subsequent description of the preconditioner . These 
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matrices are based on the bases defined above. They are 
expressed by the following relations: 

M1„ = /^,(A) -^(A)dS A (32) 

Aer 

M2 (J = J Vj (A)div A (<p j )dS A (33) 
Atr 

M3 ij= j(^(A)/sz(A)) -9,(A)dS A (34) 
Asr 

where z(A) is a vector with unit norm, normal to the 
surface T at the point A and oriented outward. 

M4 U = J M / j (A)9 l (A)dS A (35) 

Aer 

M5 g = Je J (A)9 i (A)dS A (36) 

Aer 

M6 (j =(^t(ej)J v (37) 



The steps which make it possible to compute the 
preconditioned residual from the residual in relations 
(18) and (22) is now described. These relations use the 
preconditioner Z which is now described. In the 
description which follows, we give a decomposition of 
this matrix Z into a product of matrices. The matrices 
involved in this product are sparse matrices or inverse 
sparse matrices. Thus, according to an advantageous 
variant of the invention, use will preferably be made 
of sparse matrices rather than the matrix Z directly, 
thereby making it possible in particular to reduce the 
memory used and the computation time. 



We now denote by R the residual and S the 
25 preconditioned residual. R corresponds respectively to 
R[0] and to R[n+1] in relations (18) and (22) . S 
corresponds respectively to S[0] and to S[n+1] in 
relations (18) and (22). 
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A first step consists in using a mixed representation. 
This step will be better understood with the aid of the 
vector notations which will then be transposed into 
matrix notations . The residual R corresponds to a 
bilinear form over the space Dh defined above. This 
bilinear form is denoted p . We represent p by a vector 
field rl in Dh of zero divergence and a function ql in 
Ch of zero integral. The field rl and the function ql 
are defined by the following relation for any vector 
field v in Dh: 



The vector field rl having zero divergence, we can 
write for any function f in Ch and of zero integral: 



Subsequently, so as not to needlessly complicate the 
presentation, we shall not repeat that the scalar 
fields considered are of zero integral. 

Relations (38) and (39) are transposed in matrix 
fashion as follows : 



In relation (40) , we have used a blockwise matrix 
notation in which: 

• fc M2 is the matrix transposed from M2; 

• 0 is a zero block (vector or matrix) ; 

• Rl is the matrix representation of fl ; 

• Ql is the matrix representation of ql . 



JM(A) - v(A) dS A + Jql(A) div A (v) dS A = p(v) 



(38) 




(39) 




(40) 
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We recall that it is not necessary to invert the 
blockwise defined matrix in order to determine Rl and 
Ql from R. Specifically, this matrix is a sparse 
matrix, that is to say one which contains many zero 
terms. The procedure for solving such a system is well 
known to a person skilled in the art. It is recalled in 
"Handbook of Numerical Analysis Vol. II, Mixed and 
hybrid methods" pp 523-640 by J.E. Roberts et J.-M. 
Thomas - Ed. J.L. Lions, North-Holland, 1991. This 
document forms an integral part of the description. 

A second step consists in projecting fl and ql . These 
projections are manifested by the following relations 
in vector notation: 

r2 = £>oh(r^z) (41) 
q2 = #? Sh (ql) (42) 



where : 

• is the projection operator in Dh; 

• ^ sh is the projection operator in Sh; 

• z is a unit vector, normal to T and oriented 
outward . 

These projections defined by relations (41) and (42) 
are manifested in matrix fashion by: 



(m 0 YR2) _ (M3 

[o M5AQ2j~lM4QlJ (43) 



The solution procedure for relation (43) is similar to 
the solution procedure for relation (40) insofar as the 
matrices Ml and M5 are sparse matrices. It may be noted 
that this solution procedure is even easier than that 
for relation (40) insofar as the matrices Ml and M5 are 
moreover symmetric, positive definite and well- 
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conditioned. We thus determine R2 and Q2 from Rl and 
Ql. 

A third step consists in combining r2 and q2 using the 
5 following relation: 

F3 = r2 -c"^7i(q2) (44) 

This third step is manifested in matrix fashion by: 

10 

R3 = R2-M6Q2 (45) 

We call J the vector mapping which results from the 

composition of the three steps described above. This 
15 mapping is defined by: 

j(r) = r3 (46) 

This relation (46) transposes in matrix notations to: 

20 

J R = R3 (47) 

where J denotes the matrix corresponding to the vector 
mapping J , the matrix J being defined by the following 
25 product: 

, ft m K / M1 0VVM3 0YM1 M2YY^ IM 
J = (1 - M Ho MSJ U M4A-M2 0 J [oj (48) 

It is not of course necessary to compute J directly, 
30 since the decomposition (48) into a product of sparse 
matrices and of inverses of sparse matrices is simpler 
to use. Stated otherwise, the matrix J is computed 
implicitly during its use. 
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We denote by Z the vector mapping which corresponds to 
the matrix of the preconditioner Z. The precondi tioner 
Z is defined via the vector mapping J by: 

z = j' o m o j (49) 

5 

where 

• j* is the adjoint of j; 

• o is the composition operator for composing mappings ; 
10 • m is the linear mapping Dh— > Dh* associated with the 

bilinear form m. 

Relation (42) is manifested in matrix fashion by: 
Z=' J M J (50) 

15 

Reference is now made to figures 7 and 8 in which the 
performance of an algorithm with a preconditioner 
according to the invention is illustrated in comparison 
20 with the known techniques not using a preconditioner. 

In figure 7, the curve 7 0 represents the function 

U s t 0 J|J i n the absence of a preconditioner. The 

curve 71 represents the function Vfl y in the 

25 presence of the preconditioner described above. 

In figure 8, the curve 80 represents the function 

in the absence of a preconditioner. The 

curve 81 represents the function in the 

30 presence of the preconditioner described above. 



It is found that if a conventional convergence 
criterion T) = 10 -4 is taken in relation (25), 50 
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iterations are sufficient with the preconditioner . With 
a nonpreconditioned algorithm, a sufficient accuracy is 
not reached in 200 iterations. 

5 These numerical simulations have made it possible to 
show that the use of a preconditioner according to the 
invention makes it possible to speed up and to 
stabilize the iterative algorithms. 

10 Of course, the invention is not limited to the example 
used to describe it. It is possible in particular to 
use other basis functions or a different mesh from 
those taken by way of example. 

15 The invention generalizes to any other discretization 
of the space Dh. If we take a space other than the 
Raviart-Thomas space for Dh, then the spaces Ch and Sh 
are replaced respectively by: 



Stated otherwise, 

• Ch is a minimal finite element space such that the 
divergence of the elements of Dh lies in Ch; 

2 5 • Sh is a maximal finite element space such that the 

curl of the elements of Sh lies in Dh. 

A main application of the invention is found in antenna 
design tools, but the invention is not limited to this 
30 application alone. The invention applies also of course 
to any simulation tool based on the computation of the 
field radiated by a conductor. Mention may be made in 
particular of the computation of radar cross sections 
(RCS) of objects whose geometrical properties are 

3 5 known . 




(51) 



20 




(52) 



It should be also be noted that the preconditioning 
technique described in cases of an iterative method 
also applies to other fast numerical methods. These 
fast methods are based on an iterative solution 
5 procedure, but only the terms which are useful to the 
matrix- vector products are computed. Thus, of the order 
of N x log(N) elements of the impedance matrix are 
computed, instead of N 2 elements according to the 
conventional iterative techniques. Stated otherwise, 

10 the impedance matrix is computed implicitly. The use of 
the precondi tioner according to the invention in these 
fast methods is achieved without difficulty. These 
methods are beneficial in respect of objects of large 
size, that is to say for N large. Such is the case in 

15 particular for so-called on-structure antenna 
simulations. Mention may be made in particular of the 
Multilevel Multipole Methods (or Fast Multilevel 
Multipole Methods) (FMM) and the Adaptive Integral 
Methods (AIM) . A person skilled in the art will be able 

2 0 to investigate technical elements regarding: 

• fast methods in general in "Fast Solution Methods In 
Electromagnetics" by W.C. Chew, J.-M. Jin, C.-C. Lu, 
E. Michielssen, J.M. Song - IEEE Trans, on Antennas 

25 and Propagation, 45 ( 3 ): 533-543 , March 1997; 

• Multilevel Multipole Methods in "Multilevel Fast 
Multipole Algorithm For Electromagnetic Scattering By 
Large Complex Objects" by J.M. Song, C.-C. Lu, W.C. 
Chew, S.W. Lee - IEEE Trans on Antennas and 

30 Propagation, 45 ( 10 ): 1488-1493 , October 1997; 

• Adaptive Integral Methods in "AIM: Adaptative 
Integral Method for Solving Large Scale 
Electromagnetic Scattering And Radiation Problems" by 
E. Bleszynski, M. Bleszynski, T. Jaroszewicz - Radio 

35 Science, 31 (5) : 1225-1251, 1996; 

insofar as these documents form an integral part of the 
description. 
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The invention extends without difficulty to objects 
(antennas or targets) comprising dielectric materials. 
In this case, equivalent electric and magnetic currents 
are sought on each interface. The interaction matrix 
5 between these currents comprises diagonal blocks of the 
same type as the matrix M described above. A 
preconditioner for the interaction matrix is therefore 
obtained by considering the blockwise diagonal matrix, 
whose blocks are of the type of the matrix Z described 
1 0 above . 



